!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!> One-dimensional spectral basis and associated tensor product bases
!!
!! \n
!!
!! This module provides a representation of a one-dimensional nodal
!! basis associated with a set of Gauss-Lobatto points. Tensor product
!! multidimensional extensions are also provided. Notice that the
!! following assumptions hold:
!! <ul>
!!  <li> the basis is Lagrangian, so that \f$l_i(x_j)=\delta_{ij}\f$
!!  <li> since we use Gauss-Lobatto points, the values of the basis
!!  functions at the element boundaries are either zero or one, as a
!!  special case of the Lagrangian property
!!  <li> the only nontrivial aspect of the basis functions is in fact
!!  the value of the corresponding derivatives.
!! </ul>
!! Nevertheless, when dealing with functions of different order, one
!! might need the values at points different from the Lagrangian
!! nodes: such values can be evaluated with polynomial interpolation,
!! using the Gauss-Lobatto points as interpolation nodes, as done in
!! the \c lagr_base subroutine.
!!
!! Concerning the polynomial spaces considered here, we include
!! information for both the standard piecewise polynomial basis and
!! for the Raviart-Thomas finite element space, each base with the
!! relative Gauss-Lobatto points and weights.
!<----------------------------------------------------------------------
module mod_tps_base

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_octave_io, only: &
   mod_octave_io_initialized, &
   write_octave,   &
   read_octave,    &
   read_octave_al, &
   locate_var
 
 use mod_numquad, only: &
   mod_numquad_initialized, &
   get1d_gausslobatto_nodes

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_tps_base_constructor, &
   mod_tps_base_destructor,  &
   mod_tps_base_initialized, &
   t_tps_base, new_tps_base, clear, &
   write_octave

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members

 !> Extended multiindex: this type is used to separate the extended
 !! and the nonextended components of a multiindex
 !! \f$w_{\underline{i}^{k+}}\f$. See \c t_tps_base for further details.
 type t_xmld
  integer :: xidx
  integer, allocatable :: mldx(:)
  integer, allocatable :: mldxp(:)
 end type t_xmld

 !> Base
 !!
 !! The base is essentially one-dimensional; however, it also includes
 !! index arrays to deal with multidimensional aspects of the problem,
 !! specifically multiindexes. It also includes some constants that
 !! are often required, as well as multidimensional quadrature
 !! weights.
 !!
 !! To allow the use of Raviart-Thomas basis functions, the type \c
 !! t_tps_base in fact includes <em>two</em> one-dimensional basis
 !! functions with degrees \c k and \c kx (extended basis),
 !! respectively, and the relative tensor product data. Notice that we
 !! don't assume here any relation between these two polynomial
 !! degrees, although in practice one typically has \f$k_x=k+1\f$ (for
 !! Raviart-Thomas) and \f$k_x=k\f$ (for DG).
 type t_tps_base

  integer :: d  !< spatial dimension
  integer :: pk !< number of basis functions
  integer :: k  !< polynomial degree, <tt>pk-1</tt>
  real(wp), allocatable :: xgl(:) !< Gauss-Lobatto nodes
  real(wp), allocatable :: wgl(:) !< Gauss-Lobatto weights
  real(wp), allocatable :: dphi(:,:) !< derivatives at the GL nodes
  real(wp), allocatable :: wdphi(:,:) !< \f$w_j\phi'(x^G_j)\f$

  !> Dimensions of the \f$d\f$-dimensional basis pk**d
  integer :: pkd
  !> Number of points per side: pk**(d-1)
  integer :: spkd
  !> On each row, the multiindex associated with a scalar space degree
  !! of freedom
  integer, allocatable :: mldx(:,:)
  !> Boundary degrees of freedom ( 2 , d , spkd )
  !!
  !! List the boundary degrees of freedom of each element. The first
  !! index selects "start" and "end" sides, the second index selects
  !! the direction, the third index lists the degrees of freedom.
  !!
  !! Notice that the degrees of freedom on the start and end sides are
  !! listed with a consistent ordering, so that the end degree of
  !! freedom of one element is collocated with the start degree of
  !! freedom the neighbouring one.
  integer, allocatable :: s_dofs(:,:,:)
  !> Multidimensional quadrature weigths (pkd)
  real(wp), allocatable :: wgld(:)
  !> Multidimensional side quadrature weigths ( spkd , d )
  real(wp), allocatable :: swgld(:,:)

  !> Extended polynomial order
  integer :: kx
  integer :: pkx !< number of extended basis functions
  real(wp), allocatable :: xglx(:) !< extended Gauss-Lobatto nodes
  real(wp), allocatable :: wglx(:) !< extended Gauss-Lobatto weights
  !> mixed interpolation matrix
  real(wp), allocatable :: phix(:,:)
  !> mixed derivative matrix: scalar basis evaluated at the extended
  !! GL nodes
  real(wp), allocatable :: dphix(:,:)
  !> mixed derivative matrix including the \c wglx weights
  real(wp), allocatable :: wdphix(:,:)

  !> Dimensions of the \f$d\f$-dimensional basis pkx*pk**(d-1)
  integer :: pkdx
  !> This field is analogous to \c mldx, but it refers to extended
  !! indexes.
  !!
  !! Given the extended multiindex \f$w_{\underline{i}^{k+}}\f$, the
  !! multiindex component \c mldx is given by
  !! \f$w_{\underline{i}^{k+}_{\backslash k}}\f$, while the extended
  !! component is the scalar \f$w_{\underline{i}^{k+}(k)}\f$. The
  !! complete multiindex is stored in \c mldxp
  type(t_xmld), allocatable :: mldxx(:,:)
  !> This field is analogous to \c s_dofs, but takes into account the
  !! extended basis.
  !!
  !! Notice that we are only interested here in <em>normal
  !! components</em>, so for each direction we have to check only the
  !! corresponding component, not all the \c d components. As a
  !! result, the size of \c s_dofsx is the same as the size of \c
  !! s_dofs (although the values are different).
  integer, allocatable :: s_dofsx(:,:,:)
  !> Multidimensional extended quadrature weights (pkdx).
  !!
  !! Notice that we need two indexes to store the weights
  !! \f$w_{\underline{i}^{k+}}\f$: the first index is the scalar index
  !! corresponding to \f$\underline{i}\f$, ranging from 1 to \c pkdx
  !! (this index is analogous to the index used in \c wgld), while the
  !! second index is \f$k\f$. Even though such a representation is
  !! somewhat redundant, because the weights are simply permuted by a
  !! change in \f$k\f$, it simplifies access to the correct weights.
  real(wp), allocatable :: wgldx(:,:)

  !> Shifts used by dotmap
  integer, allocatable, private :: s(:), sx(:,:)
  integer, allocatable, private :: ii(:) !< work array used by dotmap
 contains
  !> This map provides the start, stride and end indexes for slicing a
  !! local dof array.
  !!
  !! Given a multiindex \f$[i_1\,,\ldots,\,i_d]\f$ for an array with
  !! dimensions \f$[d_1\,,\ldots,\,d_d]\f$, indexing works as
  !! \f{displaymath}{
  !!  i_1 + (i_2-1)s_2 + \ldots + (i_d-1)s_d
  !! \f}
  !! where the shifts are defined as
  !! \f{displaymath}{
  !!  s_k = \prod_{h=1}^{k-1} d_h.
  !! \f}
  !! Now consider the multiindex \f$[i_1\,,\ldots,\,i_{j-1}\,, k\,,
  !! i_{j+1}\,,\ldots,\,i_d]\f$: this corresponds to
  !! \f{displaymath}{
  !!  i_1 + (i_2-1)s_2 + \ldots + (i_{j-1}-1)s_{j-1} +
  !!  (k-1)s_j + (i_{j+1}-1)s_{j+1} + \ldots + (i_d-1)s_d.
  !! \f}
  !! If we want to access the element corresponding to
  !! \f$k=k_1:k_2\f$, this corresponds to
  !! \f{displaymath}{
  !!  \begin{array}{rl}
  !!   {\rm start} &  
  !!  i_1 + (i_2-1)s_2 + \ldots + (i_{j-1}-1)s_{j-1} +
  !!  (k_1-1)s_j + (i_{j+1}-1)s_{j+1} + \ldots + (i_d-1)s_d \\
  !!   {\rm end} & 
  !!  i_1 + (i_2-1)s_2 + \ldots + (i_{j-1}-1)s_{j-1} +
  !!  (k_2-1)s_j + (i_{j+1}-1)s_{j+1} + \ldots + (i_d-1)s_d \\
  !!   {\rm stride} & s_j
  !!  \end{array}
  !! \f}
  !!
  !! This map return, \f$j=1,\ldots,d\f$, the start, end and stride
  !! indexes for a whole slice, i.e. \f$k_1=1\f$ and \f$k_2=p_k\f$.
  procedure, pass(b) :: dotmap
  !> Same as \c dotmap but for the extended basis.
  !!
  !! Notice however that, since the number of degrees of freedom is
  !! not the same in all directions for the extended basis, the
  !! indexes returned by this function can not be used for an
  !! arbitrary component of a vector field expressed with respect to
  !! the extended basis itself. Indeed, such indexes can be used only
  !! to access the \f$j\f$-th component, being \c j the last argument
  !! passed to \c dotmapx.
  procedure, pass(b) :: dotmapx
 end type t_tps_base

 ! private members

! Module variables

 ! public members
 logical, protected ::               &
   mod_tps_base_initialized = .false.
 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_tps_base'

 interface write_octave
   module procedure write_base, write_xmld
 end interface write_octave

 interface clear
   module procedure clear_base
 end interface

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_tps_base_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
       (mod_kinds_initialized.eqv..false.)    .or. &
       (mod_numquad_initialized.eqv..false.)  ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_tps_base_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   mod_tps_base_initialized = .true.
 end subroutine mod_tps_base_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_tps_base_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_tps_base_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_tps_base_initialized = .false.
 end subroutine mod_tps_base_destructor

!-----------------------------------------------------------------------
 
 !> New Gauss-Lobatto Lagrangian basis, degree k
 pure subroutine new_tps_base(b,k,d,kx)
  integer, intent(in) :: k, d
  integer, intent(in), optional :: kx !< using \c k if absent
  type(t_tps_base), intent(out) :: b
 
  integer :: i, j, h, ii, jj, counter(d)
  real(wp), allocatable :: xgl2(:,:) ! used to match the interface
  character(len=*), parameter :: &
    this_sub_name = 'new_tps_base'

   !--------------------------------------------------------------------
   ! 1) One-dimensional basis
   b%k  = k
   b%pk = k+1
   call get1d_gausslobatto_nodes(2*b%pk-3,xgl2,b%wgl)
   allocate(b%xgl(b%pk)); b%xgl = xgl2(1,:)
   deallocate(xgl2)

   allocate(b%dphi(b%pk,b%pk))
   call lagr_base(b%xgl,b%xgl, dphi=b%dphi)
   allocate(b%wdphi(b%pk,b%pk))
   do i=1,b%pk
     b%wdphi(i,:) = b%wgl*b%dphi(i,:)
   enddo

   !--------------------------------------------------------------------
   ! 2) Multidimensional basis
   b%d = d
   b%pkd  = b%pk**b%d
   b%spkd = b%pk**(b%d-1)
   allocate(b%s(b%d)); b%s = b%pk**(/(i-1, i=1,b%d)/)
   allocate(b%mldx(b%pkd,b%d))
   do i=1,b%pkd
     ii = i-1 ! it's easy with zero based indexes
     do j=b%d,1,-1
       b%mldx(i,j) = ii/b%s(j) ! using integer division
       ii = ii-b%mldx(i,j)*b%s(j)
     enddo
   enddo
   b%mldx = b%mldx + 1

   allocate(b%s_dofs(2,b%d,b%spkd))
   counter = 0
   do i=1,b%pkd
     do j=1,b%d
       if(b%mldx(i,j).eq.1) then ! found a start side dof, dir. j
         counter(j) = counter(j)+1
         b%s_dofs(1,j,counter(j)) = i
         b%s_dofs(2,j,counter(j)) = i+b%s(j)*(b%pk-1)
       endif
     enddo
   enddo

   allocate(b%wgld(b%pkd),b%swgld(b%spkd,b%d))
   do i=1,b%pkd
     b%wgld(i) = product(b%wgl( b%mldx(i,:) ))
   enddo
   do j=1,b%d
     do i=1,b%spkd
       ! get the corresponding multiindex (for instance, from ie1)
       counter = b%mldx(b%s_dofs(2,j,i),:) ! used as temporary
       b%swgld(i,j) = product(b%wgl(counter))/b%wgl(counter(j))
     enddo
   enddo

   !--------------------------------------------------------------------
   ! 3) Extended one-dimensional basis
   if(present(kx)) then
     b%kx = kx
   else
     b%kx = k
   endif
   b%pkx = b%kx+1
   call get1d_gausslobatto_nodes(2*b%pkx-3,xgl2,b%wglx)
   allocate(b%xglx(b%pkx)); b%xglx = xgl2(1,:)
   deallocate(xgl2)

   allocate(b%phix(b%pk,b%pkx),b%dphix(b%pk,b%pkx))
   call lagr_base(b%xglx,b%xgl , phi=b%phix,dphi=b%dphix)
   allocate(b%wdphix(b%pk,b%pkx))
   do i=1,b%pk
     b%wdphix(i,:) = b%wglx*b%dphix(i,:)
   enddo

   !--------------------------------------------------------------------
   ! 4) Extended multidimensional basis
   b%pkdx = b%pkx*(b%pk**(b%d-1))
   allocate(b%sx(b%d,b%d))
   do h=1,b%d
     do i=1,h
       b%sx(i,h) = b%pk**(i-1)
     enddo
     do i=h+1,b%d
       b%sx(i,h) = b%pk**(i-2) * b%pkx
     enddo
   enddo
   allocate(b%mldxx(b%pkdx,b%d))
   do h=1,b%d
     do i=1,b%pkdx
       allocate(b%mldxx(i,h)%mldx(b%d-1))
       allocate(b%mldxx(i,h)%mldxp(b%d))
       ii = i-1 ! it's easy with zero based indexes
       do j=b%d,1,-1
         jj = ii/b%sx(j,h) ! using integer division
         ii = ii-jj*b%sx(j,h)
         ! Compared to the mldx case, we have to differentiate j=h
         if(j.gt.h) then
           b%mldxx(i,h)%mldx(j-1) = jj+1 ! restore one based indexes
           b%mldxx(i,h)%mldxp(j) = b%mldxx(i,h)%mldx(j-1)
         elseif(j.eq.h) then
           b%mldxx(i,h)%xidx      = jj+1
           b%mldxx(i,h)%mldxp(j) = b%mldxx(i,h)%xidx
         else
           b%mldxx(i,h)%mldx( j ) = jj+1
           b%mldxx(i,h)%mldxp(j) = b%mldxx(i,h)%mldx( j )
         endif
       enddo
     enddo
   enddo

   allocate(b%wgldx(b%pkdx,b%d))
   do j=1,b%d
     do i=1,b%pkdx
       b%wgldx(i,j) = product(b%wgl( b%mldxx(i,j)%mldx )) &
                           * b%wglx( b%mldxx(i,j)%xidx )
     enddo
   enddo

   allocate(b%s_dofsx(2,b%d,b%spkd))
   counter = 0
   ! notice: whether a given extended multiindex is on the boundary or
   ! not depends also on the considered direction; for each direction,
   ! anyway, we are only interested in degrees of freedom of the
   ! corresponding vector components, as long as we only need normal
   ! fluxes.
   do h=1,b%d
     do i=1,b%pkdx
       if(b%mldxx(i,h)%mldxp(h).eq.1) then ! see s_dofs
         counter(h) = counter(h)+1
         b%s_dofsx(1,h,counter(h)) = i
         ! In the next row we exploit the fact that along the h-th
         ! direction and for the h-th component the number of degrees
         ! of freedom is pkx (this would be different along different
         ! directions).
         b%s_dofsx(2,h,counter(h)) = i+b%sx(h,h)*(b%pkx-1)
       endif
     enddo
   enddo

   ! Other fields
   allocate(b%ii(b%d))

 end subroutine new_tps_base
 
!-----------------------------------------------------------------------

 pure subroutine clear_base(base)
  type(t_tps_base), intent(out) :: base
  character(len=*), parameter :: &
    this_sub_name = 'clear_base'
   
   ! nothing to do
 end subroutine clear_base
 
!-----------------------------------------------------------------------
 
 pure subroutine lagr_base(x,xn, phi,dphi,ddphi)
  real(wp), intent(in) :: xn(:) !< nodes
  real(wp), intent(in) ::  x(:) !< coordinates
  real(wp), intent(out), optional :: phi(:,:), dphi(:,:), ddphi(:,:)
 
  integer :: i, j, k, l, pk
  real(wp) :: prod(size(x)), sum(size(x))
  character(len=*), parameter :: &
    this_sub_name = 'lagr_base'

  pk = size(xn)

  ! phi
  if(present(phi)) then
    phi = 1.0_wp
    do i=1,pk ! basis functions
      do j=1,pk ! nodes
        if(j.ne.i) & ! Lagrange polynomials
          phi(i,:) = phi(i,:) * (x-xn(j))/(xn(i)-xn(j))
      enddo
    enddo
  endif

  ! dphi
  if(present(dphi)) then
    dphi = 0.0_wp
    do i=1,pk ! basis functions
      do k=1,pk ! nodes
        if(k.ne.i) then
          prod = 1.0_wp
          do j=1,pk ! nodes
            if((j.ne.i).and.(j.ne.k)) &
              prod = prod * (x-xn(j))/(xn(i)-xn(j))
          enddo
          dphi(i,:) = dphi(i,:) + prod/(xn(i)-xn(k))
        endif
      enddo
    enddo
  endif
 
  ! ddphi
  if(present(ddphi)) then
    ddphi = 0.0_wp
    do i=1,pk ! basis functions
      do k=1,pk ! nodes
        if(k.ne.i) then
          sum = 0.0_wp
          do l=1,pk
            if((l.ne.i).and.(l.ne.k)) then
              prod = 1.0_wp
              do j=1,pk ! nodes
                if((j.ne.i).and.(j.ne.k).and.(j.ne.l)) &
                  prod = prod * (x-xn(j))/(xn(i)-xn(j))
              enddo
              sum = sum + prod/(xn(i)-xn(l))
            endif
          enddo
          ddphi(i,:) = ddphi(i,:) + sum/(xn(i)-xn(k))
        endif
      enddo
    enddo
  endif
 
 end subroutine lagr_base
 
!-----------------------------------------------------------------------
 
 !> Argument \c b should be intent(in), except that we use the work
 !! array \c ii.
 !!
 !! \todo This subroutine seems to be performance critical; any
 !! improvement would be useful.
 !!
 !! When using OpenMP, the work array ii would be shared by all the
 !! threads, leading to race conditions. In such a case, the work
 !! array must be passed by the caller.
 pure subroutine dotmap(b,start,end,stride,i,j &
 !$                     , ii &
                       )
  class(t_tps_base), intent(inout) :: b
  integer, intent(in) :: i(:), j
  !$ integer, intent(out) :: ii(:)
  integer, intent(out) :: start, end, stride
 
  character(len=*), parameter :: &
    this_sub_name = 'dotmap'

   !$ openmp_if: if(.false.) then
   b%ii = i-1
   b%ii(j) = 0
   start  = dot_product(b%ii,b%s) + 1
   !$ else ! OpenMP version
   !$   ii    = i-1
   !$   ii(j) = 0
   !$   start = dot_product(ii,b%s) + 1
   !$ endif openmp_if
   end    = start + (b%pk-1)*b%s(j)
   stride = b%s(j)

 end subroutine dotmap
 
!-----------------------------------------------------------------------

 !> See \c dotmap
 pure subroutine dotmapx(b,start,end,stride,i,j &
 !$                      , ii &
                        )
  class(t_tps_base), intent(inout) :: b
  integer, intent(in) :: i(:), j
  !$ integer, intent(out) :: ii(:)
  integer, intent(out) :: start, end, stride
 
  character(len=*), parameter :: &
    this_sub_name = 'dotmapx'

   !$ openmp_if: if(.false.) then
   b%ii = i-1
   b%ii(j) = 0
   start  = dot_product(b%ii,b%sx(:,j)) + 1
   !$ else ! OpenMP version
   !$   ii    = i-1
   !$   ii(j) = 0
   !$   start = dot_product(ii,b%sx(:,j)) + 1
   !$ endif openmp_if
   end    = start + (b%pkx-1)*b%sx(j,j)
   stride = b%sx(j,j)

 end subroutine dotmapx
 
!-----------------------------------------------------------------------

 subroutine write_xmld(xmld,var_name,fu)
  integer, intent(in) :: fu
  type(t_xmld), intent(in) :: xmld(:,:)
  character(len=*), intent(in) :: var_name
 
  integer :: i, j
  character(len=*), parameter :: &
    this_sub_name = 'write_xmld'

   write(fu,'(a,a)')    '# name: ',var_name
   write(fu,'(a)')      '# type: struct'
   write(fu,'(a)')      '# length: 3' ! number of fields

   ! field 01 : xidx
   write(fu,'(a)')      '# name: xidx'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a,i7)')   '# rows: ',size(xmld,1)
   write(fu,'(a,i7)')   '# columns: ',size(xmld,2)
   do j=1,size(xmld,2)
     do i=1,size(xmld,1)
       call write_octave(xmld(i,j)%xidx,'<cell-element>',fu)
     enddo
   enddo

   ! field 02 : mldx
   write(fu,'(a)')      '# name: mldx'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a,i7)')   '# rows: ',size(xmld,1)
   write(fu,'(a,i7)')   '# columns: ',size(xmld,2)
   do j=1,size(xmld,2)
     do i=1,size(xmld,1)
       call write_octave(xmld(i,j)%mldx,'r','<cell-element>',fu)
     enddo
   enddo

   ! field 03 : mldxp
   write(fu,'(a)')      '# name: mldxp'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a,i7)')   '# rows: ',size(xmld,1)
   write(fu,'(a,i7)')   '# columns: ',size(xmld,2)
   do j=1,size(xmld,2)
     do i=1,size(xmld,1)
       call write_octave(xmld(i,j)%mldxp,'r','<cell-element>',fu)
     enddo
   enddo
 
 end subroutine write_xmld

!-----------------------------------------------------------------------

 subroutine write_base(b,var_name,fu)
  integer, intent(in) :: fu
  type(t_tps_base), intent(in) :: b
  character(len=*), intent(in) :: var_name
 
  character(len=*), parameter :: &
    this_sub_name = 'write_base'

   write(fu,'(a,a)')    '# name: ',var_name
   write(fu,'(a)')      '# type: struct'
   write(fu,'(a)')      '# length: 26' ! number of fields

   ! field 01 : pk
   write(fu,'(a)')      '# name: pk'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%pk,'<cell-element>',fu)

   ! field 02 : k
   write(fu,'(a)')      '# name: k'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%k,'<cell-element>',fu)

   ! field 03 : xgl
   write(fu,'(a)')      '# name: xgl'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%xgl,'r','<cell-element>',fu)

   ! field 04 : wgl
   write(fu,'(a)')      '# name: wgl'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%wgl,'r','<cell-element>',fu)

   ! field 05 : dphi
   write(fu,'(a)')      '# name: dphi'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%dphi,'<cell-element>',fu)

   ! field 06 : wdphi
   write(fu,'(a)')      '# name: wdphi'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%wdphi,'<cell-element>',fu)

   ! field 07 : d
   write(fu,'(a)')      '# name: d'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%d,'<cell-element>',fu)

   ! field 08 : pkd
   write(fu,'(a)')      '# name: pkd'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%pkd,'<cell-element>',fu)

   ! field 09 : spkd
   write(fu,'(a)')      '# name: spkd'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%spkd,'<cell-element>',fu)

   ! field 10 : mldx
   write(fu,'(a)')      '# name: mldx'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%mldx,'<cell-element>',fu)

   ! field 11 : s_dofs
   write(fu,'(a)')      '# name: s_dofs'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%s_dofs,'<cell-element>',fu)

   ! field 12 : wgld
   write(fu,'(a)')      '# name: wgld'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%wgld,'r','<cell-element>',fu)

   ! field 13 : swgld
   write(fu,'(a)')      '# name: swgld'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%swgld,'<cell-element>',fu)

   ! field 14 : kx
   write(fu,'(a)')      '# name: kx'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%kx,'<cell-element>',fu)

   ! field 15 : pkx
   write(fu,'(a)')      '# name: pkx'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%pkx,'<cell-element>',fu)

   ! field 16 : xglx
   write(fu,'(a)')      '# name: xglx'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%xglx,'r','<cell-element>',fu)

   ! field 17 : wglx
   write(fu,'(a)')      '# name: wglx'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%wglx,'r','<cell-element>',fu)

   ! field 18 : phix
   write(fu,'(a)')      '# name: phix'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%phix,'<cell-element>',fu)

   ! field 19 : dphix
   write(fu,'(a)')      '# name: dphix'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%dphix,'<cell-element>',fu)

   ! field 20 : wdphix
   write(fu,'(a)')      '# name: wdphix'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%wdphix,'<cell-element>',fu)

   ! field 21 : pkdx
   write(fu,'(a)')      '# name: pkdx'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%pkdx,'<cell-element>',fu)

   ! field 22 : mldxx
   write(fu,'(a)')      '# name: mldxx'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%mldxx,'<cell-element>',fu)

   ! field 23 : wgldx
   write(fu,'(a)')      '# name: wgldx'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%wgldx,'<cell-element>',fu)

   ! field 24 : s_dofsx
   write(fu,'(a)')      '# name: s_dofsx'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%s_dofsx,'<cell-element>',fu)

   ! field 25 : s
   write(fu,'(a)')      '# name: s'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%s,'r','<cell-element>',fu)

   ! field 26 : sx
   write(fu,'(a)')      '# name: sx'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(b%sx,'<cell-element>',fu)

 end subroutine write_base

!-----------------------------------------------------------------------

end module mod_tps_base

